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ABSTRACT 

The model for pulsar wind nebulae (PWNe) as the result of the magnetohydrodynamic (MHD) 
downstream flow from a shocked, relativistic pulsar wind has been successful in reproducing many 
features of the nebulae observed close to the central pulsars. However, observations of well-studied 
young nebulae like the Crab Nebula, 3C 58, and G21.5-0.9 do not show the toroidal magnetic field on a 
larger scale that might be expected in the MHD flow model; in addition, the radial variation of spectral 
index due to synchrotron losses is smoother than expected in the MHD flow model. We find that 
pure diffusion models can reproduce the basic data on nebular size and spectral index variation for the 
Crab, 3C 58, and G21.5-0.9. Most of our models use an energy independent diffusion coefficient; power 
law variations of the coefficient with energy are degenerate with variation in the input particle energy 
distribution index in the steady state, transmitting boundary case. Energy dependent diffusion is a 
possible reason for the smaller diffusion coefficient inferred for the Crab. Monte Carlo simulations of 
the particle transport allowing for advection and diffusion of particles suggest that diffusion dominates 
over much of the total nebular volume of the Crab. Advection dominates close to the pulsar and is 
likely to play a role in the X-ray half-light radius. The source of diffusion and mixing of particles is 
uncertain, but may be related to the Rayleigh- Taylor instability at the outer boundary of a young 
PWN or to instabilities in the toroidal magnetic field structure. 
Subject headings: ISM: individual (Crab Nebula) 
stars: winds, outflows 



-ISM: supernova remnants — pulsars: general 



1. INTRODUCTION 

The finding of the diminishing size of the Crab syn- 
chrotron nebula with increasing frequency supports the 
picture that energetic electrons are injected in the vicin- 
ity of the pulsar and l ose energy to synchrotron radiatio n 
in the larger nebula (|Wilsonlll97l iRees fe Gunnlll974f ). 
The radial emission profiles were first modeled as a cen- 
tral particle source with diff u sion into the larger n ebula 
(|Grattonlll97ilWilsonlll97l . IRees fc Guniil (fl97l spec- 
ified a termination shock in the pulsar wind as the source 
of the particle acceleration and viewed the outer part of 
the nebula as a place where the pulsar magnetic field 
winds up and the flow decelerates. This view was put on 
a firmer basis bv lKennel fe Coronitil (|1984alfbL hereafter 
KC), who calculated the conditions at the relativistic 
MHD shock and followed the time independent down- 
stream flow of fields and particles. This 1-dimensional 
model with advected particles was able to reproduce the 
observed sizes of the optical and X-ray emission in the 
Crab Nebula, but did not address the radio emission. 

High resolution imaging of the Crab with the Hubble 
Telescope at optical wavelengths and Chandra at X-ray 
wavelengths shows an active system of toroidal filaments 
and jets close to the pulsar in li ne with the e xpected 
position of the termination shock ((Hester 2008]). These 
observations have motivated 2-dimensional MHD sim- 
ulations, which allow f or a polar angle dependence of 
the pulsar wind power (jKomissarov fe Lvubarskvl 120031 : 
iDel Zanna et al.l 12004) . In these models, the wind is 
stronger in the equatorial plane, producing the toroidal 
filaments. Hoop stresses in the shocked flow bring ma- 
terial back to the axis to form the jets. Current models 
for the filaments are able to reproduce many aspects of 



the filaments (see [Buccian tlriil 1201 ll for a review), in- 
cluding the integrated spect rum of the Crab from ra- 
dio to TeV dVolpi et al.ll2008ri and the time variability of 
the inner structure (iVolpi et al.ll200llCamus et al.l[200l 
iKomissarov fe Lvutikovl I2011D . The jet-torus structure 
near the pulsar has bee n commonly observed in X-ra y im- 
ages of pulsar nebulae (Kargalts ev fc Pavlovll200"8l ). and 
the structure is presumed to be a standard feature of pul- 
sar nebulae. Although models for the toroidal filaments 
are now convincing, the nature of the flow beyond the 
filaments to the edge of the nebula remains uncertain. 

The primary way of exploring the particle transport is 
to model the structure of the nebulae at different pho- 
ton energies or, equivalent ly, the structure and photon 
index distribution. The crucial point is that the par- 
ticles lose energy to synchrotron radiation as they age 
so that the photon index distribution provides a good 
test of the particle transport mechanism. Provided the 
magnetic field is not strongly varying, the spectral in- 
dex in a particular location gives information on the 
mean age of the particles. Observations of the Crab 
Nebula at optical (IVeron-Cettv fc Woltierlll99l and in- 
frared (|Temim et al.l 12006) wavelengths show a mono- 
tonic change in spectral index from the center to the 
edge of the nebula, where edge is defined by the decrease 
in surface brightness at that particular wavelength. The 
well-known bays in the Crab are asymmetric structures, 
but the spectral index at the bay edge is similar to that 
at other edges. The data do not indicate a highly asym- 
metric flow in the nebula. Thus, we assume spherical 
symmetry in our models. 

In Section 2 of this paper, we discuss issues with exist- 
ing models for the larger scale flow and, in Section 3, con- 
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sider a diffusion/ advection model for the particle propa- 
gation. The model is applied to the phase when the pul- 
sar wind nebula (hereafter PWN) is expanding into the 
freely expanding supernova gas, before a reverse shock 
front moves in due to interaction with the surroundings. 
We concentrate on comparing our models to the Crab 
Nebula, 3C 58 and G21.5-0.9 because they are the best 
observed PWNe within that phase. In Section 4, we dis- 
cuss the diffusion process. 

2. OUTER STRUCTURE OF YOUNG PULSAR NEBULAE 

Current data on the Crab Nebula convincingly show 
that there is a relativistic flow in th e interior tha t must 
slow to match the outer boundary (Hester 2008). The 
issue treated in this paper is how the particles are trans- 
ported between the inner, toroidal region and the outer 
boundary of the PWN. In the KC model, the parti- 
cles are advected with a toroidal magnetic field. Cross 
field scattering of particles is ex pected to be small (e.g., 
Ide Jager fc Diannati-At a'i 2009) , so that diffusion of par- 
ticles can be neglected. Problems with this model in 
reproducing the spectral i ndex distributions in young 
PWN e have been raised (Reynolds] 120031 : ISlane et al.l 
2004). This issue is discussed in detail in the next sec- 
tion. Here we note some other points relevant to the 
toroidal field model for the outer structure. 

In the Crab Nebula, the X-ray emission is from a re- 
gion close to the pulsar because of synchrotron burn- 
off. If one goes to optical wavelengths, where the par- 
ticles have longer lifetimes, an analysis of the polariza- 
tion shows that there are 3 — 6 magnetic elements across 
the n eb ula with possible sma ller scale structure ()Feltenl 
119741 ). iSchmidt et all <| 19791) find that magnetic struc- 
ture in th e Crab only extends down to about 20 arcsec, 
or 0.2 pc. Seward et al.l (|2006f ) presented deep Chandra 
images of the Crab, finding evidence for fingers with a 
roughly radial orientation; the spectral index structure 
implied rapid diffusion along the structures, presumably 
oriented along the magnetic field and slow diffusion per- 
pendicular to the fingers. The X-ray emitting particles 
in 3C 58 have longer lifetimes and imaging X-ray studies 
with Chandra show many magnetic loops wit hout a clear 
toroi dal structure, except close to the pulsar (|Slane et al.l 
2004); the X-ray filaments arc related to ones observed at 
radio wavelengths and with some optical filaments (ther- 
mal gas). Overall, there is little evidence for toroidal 
structure in pulsar nebulae except near the central pul- 
sars. Two PWNe that do show evidence for toroidal 
field structure are the sm all nebula around the Vela pul- 
sar (iDodson et all 120031 ) and G106.6+2.9 ([Rothes et al.l 
l2006f) : however, these objects are probably in a different 
evolutionary sta ge than the Crab, without a n unstable 
outer boundary ([Chevalier fc Reynolds! |20TTI ). 

The Crab Nebula, with ra dio spectral index j3 = 
0.299±0.009 (|Baars et al.lll977l) where flux cx V 13 , shows 
little radio spectral i ndex variation over the entire neb- 
ula, to within 0.01 (|Bietenholz et al.l [19971) . although 
PWNe themselves show a range of spectral indic es, e.g., 
3C 58 has /3 = 0.07±0.05 (jBietenholz et al.ll2001l ). There 
is no indication that the spectral index observed in the 
Crab is a universal value, so the uniformity of spectral 
index is surprising. G2 1.5-0 .9 also has a fairly uniform 
radio spectral in dex image dBietenholz fc Bartell |008) , 
as well as 3C 58 (jBietenholz et al.ll2001l) . 



These observational considerations support the view 
that, although there is clearly toroidal structure where 
the pulsar wind impacts the larger nebula, the flow is 
more radial in the outer nebula and there is evidence 
for a mixing process. There are several possible rea- 
sons for the apparent mixing of energetic particles. The 
acceleration of the supernova e.jec ta by the pulsar bub- 
ble is Rayleigh- Taylor u nstable ([Chevalierl Il977t Uunl 
Il998t iHester et all 119961 ). The structure observed in 
the thermal gas filaments in the Crab Nebula and o ther 
PWNe is likely due to this instability (jHesterl 120081) . al- 
though there are still uncertainties about how the insta- 
bility operates when the low density fluid is magn etized 
(jBucciantini et al.l 120041: iStone fc Gardinerl I2007D . As 
discussed bv lHester et al.l ([1995D . there is evidence in the 
Crab for magnetic field lines being 'draped' around ther- 
mal filaments and stretched in the radial direction. The 
Crab filament s cover the velocity range 700—1500 km s _1 
(|Hesteiil2008h . 

Another possibility is the action of instabilities oc- 
curring at or near the pulsar wind termination shoc k 
(|Begelmanlll99l ICamus et ai1l2009l; IMizuno et aT1l201lD . 
Instabilities right at the shock front may not be compat- 
ible with the apparent regular structure observed in the 
case of the newly formed nebula around the Vel a pulsar 
(jDodson et al.ll200l IChevalier fc Revnoldsll201lD . There 
may be feedback between instabilities near the termina- 
tion shock and the outer boundary. In their a xisymmet- 
ric numerical simulations, ICamus et al.l ([2009) find that 
waves and vortices in the larger nebula feed back on the 
structure at the termination shock, which in turn gener- 
ates more structure in the nebula. ICamus et al.l (2009) 
note that it is impossible to distinguish between the cause 
and the effect. In addition, some thermal gas from the 
supernova is entrained in the unstable region, explaining 
the optical filaments seen in association with nonthermal 
filaments in the Crab and 3C 58. 

These observational and theoretical considerations 
show that, although the region close to the pulsar has 
a clear toroidal structure, the larger nebula has a com- 
plex structure that includes a radi al component to the 
magnetic field. As summarized by IHesterl (|2008|) , there 
are layers of magnetic fields folded on top of each other 
such that adjoining field lines in one place move away 
from each other. Although diffusion across magnetic field 
lines is expected to be small, even a small amount of cross 
field transport could result in thorough mixing, with this 
magnetic field configuration. We thus consider models 
with radial diffusion. 

One uncertainty for the models is the degree to which 
particles are transmitted through the outer boundary 
of the PWN. The expectation is that the PWN mag- 
netic field is contained within the wind bubble, which is 
bounded by supernova ejecta in young PWNe; the outer 
boundary of the Crab Nebula shows loop structures at 
radio wavel engths that are likely to delineate the mag- 
netic field ([Hesterl [20081) . The characteristic mean free 
path for particles to cross the magnetic field is limited to 
the particle gyroradius (Bohm limit) 

where the particle energy, E e , corresponds to an 
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energy of synchrotron radiation through E e = 

(20 TeV)(£/100 fj,G)~ 1 ^ 2 E^ e y and B is the magnetic 
field. The escape time for Bohm diffusion from a region 
of size R is then (jde Jager fe Diannati-At ai 2009) 

^- 16 ' 000 (2^) 2 (l0^w) 1 

x (ioo^g) yr ' (2) 

which is long compared to the ages of the PWNe consid- 
ered here (~ 10 3 yr), particularly for particles radiating 
below X-ray wavelengths. However, there is some chance 
for escape from a narrow region close to the edge of the 
nebula. 

An observational test for the transmission of particles 
would be synchrotron emission from energetic particles 
that have left the main P WN. There has been l ittle evi- 
dence for such emission, but lBamba et al.l (|2010f ) have re- 
cently found X-ray emission around a number of PWNe, 
which they interpret as synchrotron emission from es- 
caped particles. To have X-ray emitting particles extend 
to such large radial distances requires a surprisingly low 
magnetic field strength in order to avoi d synchrotron 
losses . One of the PWNe discussed by iBamba et al.l 
(2010) is G21. 5-0.9, which is also one of the primary 
remnants treated here. However, ther e are other in- 
terpre tations for the extended emission. IBocchino et al.l 
( 2005) attributed the emission to a combination of a dust 
scattering halo plus emission from a surrounding shell; 
iMatheson fc Sa fi-Harb (201(J have shown clear evidence 
for shell emission. 

Escape of electrons and positrons from PWNe has also 
come up in the context of a possible source for fea- 
tures in the Galacti c cosmic ray spect rum of electrons 
and positrons (e.g.. iChang et al1l2008|) . However, the 
escape can be from elderly PWNe (> 10 5 yr old), and 
does not require es cape from young obj e cts like those 
discus sed here (e.g.. lMalvshev et aLll2009l) . iHinton et al.l 
()2011|) suggested the escape of particles from the Vela 
X PWN to explain the steeper particle spectrum in the 
outer parts of the nebula; however, this is an older nebula 
that has likely been affected by the reverse shock wave 
(jBlondin et al.ll200lD . Overall, particle escape in PWNe 
with ages ~ 10 3 yr does not appear likely. In Section 3.2, 
we model the effect of the outer boundary condition for 
the PWN. 

In the early diffusion models for the Crab Ne bula 
(|Grattonlll97aiWilsonl[l97llWeinberg fc SiHd[l976h . the 
diffusion coefficient was assumed to be constant with en- 
ergy, and that is the assumption that we make in most 
of our modeling. Weinberg fc Silk! (|1976l ) argued for a 
constant coefficient based on the fact that the diffusion 
length is likely to be related to the size of magnetic fil- 
aments and is much larger than the gyroradius. How- 
ever, there are reasons to consider energy dependent dif- 
fusion of the form De(E) = DqE 01 , where a is a con- 
stant. On the observational side, cosmic rays are known 
to di ffuse from the Gala xy in an energy dependent way 
(e.g.. iStrong et all 120071 and references therein). The 
data indicate a = 0.3 — 0.6, with a diffusion coefficient 
of (3 — 5) x 10 28 cm 2 s" 1 at a reference particle energy 
of 1 GeV (Str ong et al.ll2007j) . Interstellar turbulence is 



thought to play a role in the energy dependent diffusion. 
Also, there is the possibility that the diffusion coefficient 
is proportional to the particle gyroradius, as in Bohm dif- 
fusion, l eading to a = 1. This "Bohm - type " value of a is 
used b y IVan Etten fc Romanil (|2011[ ) and IHinton et al.l 
(|2011[ ) in their modeling of evolved PWNe. However, 
the diffusion length is much larger than the particle gy- 
roradius, so there is not a clear argument for a = 1 in 
the PWN case. We consider the possible effect of energy 
dependent diffusion in Section 3.1. 

3. MODELS WITH DIFFUSION 
3.1. Pure diffusion model 

Wilson (1972) first showed that the spatial and spec- 
tral distribution of the Crab Nebula in the optical could 
be explained by a diffusion model with synchrotron ra- 
diation losses. Observations of 3C 58 (Boc chino et all 
20011 ISlane et al 2004) and G21. 5-0.9 (jSlane et al.l(200a 
Safi- Harb et alll200lD taken by Chandra also gave a pho- 
ton index distribution that is similar to the optical spec- 
tral index distributio n in the Crab and is incompatibl e 
with the KC model (jRevnoldsl 12001 ISlane et al.ll200l . 
Although observations of the Crab clearly show evidence 
for a relativistic wind close to the pulsar, a pure diffu- 
sion model illustrates one limiting case of the expected 
particle trans port. W e used the pure diffusion model 
developed by iGrattonl (|1972fl to fit the spectral index 
distribution of Crab from radio to optical, and 3C 58 
and G21.5-0.9 in X-rays. We also used the model to 
calculate the half light radius of the C rab from radio to 
X-rays. This is the same model used bv lWilsonl (|1972l ) in 
his model for the Crab. The model assumes that a point 
source injects particles into an infinite space with spher- 
ical symmetry, and the injected particles follow a power 
law distribution N{E,r = 0) = KE~ P . The transport 
mechanism for the injected particles is diffusion. In or- 
der to satisfy the spherical symmetry assumption, the 
objects we discuss in this paper are young PWNe that 
are observed to have approximate spherical symmetry. 
In our initial model we have a pre-defined PWN radius 
R; when we calculate the half light radius and integrated 
spectrum, we only consider the particles that are within 
the nebular radius R, reasoning that the magnetic field 
in the freely expanding ejecta outside R is small. The 
outer nebular boundary is a transmitting boundary. The 
model also assumes that the diffusion coefficient D and 
magnetic field B inside R are constant with radius. Here 
we neglect the adiabatic expansion energy losses and as- 
sume synchrotron radiation loss is the only energy loss, so 
that dE/dt = -QE 2 , where Q = 2.37 x 1Q- 3 J3 2 er g s" 1 
in cgs units and B\ = (2/3)B 2 (|Pacholczvklll970f ). We 
examine the assumptions of point source injection, pure 
diffusion, and a transmitting boundary in Section 3.2. 
In this section we use this simple pure diffusion model 
to analyze the spectral index distribution and half light 
radius of PWNe. 

Based on these assu mptions, the n umber density dis- 
tribution N(E, r, t) is (|Grattonlll972h 

N(E,r,t) = ^E-Pf p (u,v), (3) 

where 

r 2 QE 
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and 



f p (u,v) = —= 



1 

x 



rdx. 



The lower limit v of the integral f p (u,v) is 



4-Dt 
r 2 QE 
4D 



if there is no upper limit for the injection particle energy. 
Here t is the age of the nebula. 

In order to simplify the calculation of the emission, 
we further assume that all the radiated power P of an 
electron of energy E goes into radiation of a frequency 
v corresponding to the maximum synchrotron radiation 
power. Therefore 



P{v)=C(B^) 1 ' 2 N[E( V )], 
where C is a constant, and 



E{v) 



0.29 x 3e 



1/2 



B sin 6 



1/2 



(4) 



(5) 



where m e and e are the mass and charge of an elec- 
tron, and 9 is the particle pitch angle. Here we as- 
sume Bsin6 = B± = (2/3) ^B, yielding E(v) = 
7.42 x 10- 10 (v/B 1 _) 1 / 2 erg. The spectral index distri- 
bution S(r) between frequency V\ and v 2 is given by 



Sir) 



\og[P vl (r)/P v2 {r)] 



\og(v 1 /v 2 ) 

\og{vi/v 2 ) 1/2 + \og[N tot (v 1 ,r)/N tot (v 2 ,r)} 
lag(n/i*2) 



(6) 



where N tot {v,r) is the total number of particles emit- 
ted per unit area per unit time, per unit frequency with 
frequency v and at radius r from a cent ral point source 
after integration along the line of sight. iGrattonl ()1972f ) 
assumed a point source which makes r = 0a singularity 
in Ntot {v 1 r) . We performed integrations along the line 
of sight starting at a cutoff radius, but this did not affect 
the larger scale results. 

There is a critical energy for synchrotron cooling, 
E C rit = 1/Qt, which is relevant for the number density 
distribution N(E, r, t). If E > E crit , N(E, r, t) reaches a 
steady state solution N(E,r). If E < E crit , only parti- 
cles with Ei n iuai < E/(l — QEt) contribute to the spec- 
trum at the frequency v and N{E 1 r, t) evolves with time. 
The corresponding peak synchrotron radiation frequency 
of particle with E cr u is 



Vcrit = 6 x 10 1 



10 3 yr\ / 10Q jiG 



Hz. (7) 



If the injected particles have an upper limit of energy E*, 
the lower limit v of the integral f p (u,v) changes to 



4DI 



if t > A (± - ^ 



and the critical energy becomes E*/(l + QtE*). When 
the energy range of interest satisfies E <C E* , the E* 
term in v is not important and can be neglected, which 



means we can assume the injection particle energy has 
no upper limit. A plausible estimate for E* is that the 
gyroradius of the particle equals to the termination shock 
radius, R s , or R s — R gyr o — E/eB for a relativistic elec- 
tron. We then have E* — R s eB and the corresponding 
peak synchrotron emission frequency is 



v* = 3.3 x 10 



22 



Rs 



O.lpc 



f— V 

VIOO^G J 



Hz. 



For the objects considered in this paper, electron energy 
injection with no upper limit is always a good assumption 
for frequencies below soft X-rays. 

The spectral index distribution of the system at a cer- 
tain frequency band mainly depends on the ratio rj be- 
tween diffusion distance d = (6-Dt) 1 / 2 and the nebular 
size R. At a certain frequency, the same rj = d/R = 
(6Dt/ R 2 ) 1 / 2 gives roughly the same spectral index dis- 
tribution. The p index of the injected particles would also 
affect the spectral index profile to some extent. Since we 
have good observational data for the spectral index of 
the PWNe considered here, its value can be directly de- 
termined from the observations. If the frequency band 
is in the steady state regime, then t = 1/QE, and we 
find rj oc {D / v 1 / 2 B z / 2 R 2 ) 1 / 2 . Since the nebular size 
is usually known, the spectral index distribution is de- 
termined by r\ cc {D / B 3 / 2 ) 1 / 2 in a particular frequency 
band. The diffusion coefficient D and magnetic field B 
are coupled together in the spectral index fitting; in or- 
der to get an accurate diffusion coefficient D, we need to 
know the magnetic field B of the system. One way to 
estimate B is based on the synchrotron break frequency 
in the integrated spectrum (equation [7]). However, the 
break is not a sharp feature and it is difficult to locate 
the synchrotron break frequency of th e Crab and 3C 58 
based on current o bservational data (jSlane et al.l 120081 : 
lArendt et al.ll201ll ). Another way is to model the high 
energy inverse Compton emission. The magnetic field ob- 
tained by using inverse Compton fitting is slightly differ- 
ent from the magnetic field defined in our pure diffusion 
model, because in our model we solve for the constant B 
situation which means the magnetic field B is an average 
B of the PWN over its lifetime. Our average magnetic 
field should be slightly larger than the magnetic field in- 
dicated by inverse Compton modeling. The minimum 
energy method for synchrotron emission also gives an es- 
timate for B, but it depends on an uncertain assumption. 
If the frequency band is in the non-steady state regime, 
rj = d/R = (6Dt/ 'R 2 ) 1 / 2 , where t is now the age of the 
nebula. The diffusion coefficient D and magnetic field 
B are now decoupled and rj oc (D) 1 / 2 . However, in this 
regime the particles do not suffer synchrotron losses, so 
their spectrum is not changed from the injection spec- 
trum and they do not give useful information about the 
diffusion coefficient. 

Next we consider how the spectral index profile varies 
as a function of the ratio rj. As we are mainly interested 
in the steady state case, we now consider the critical fre- 
quency vr corresponding to the case rj « (4Dt) 1 / 2 /R = 
(AD/QER 2 ) 1 / 2 = 1. It is the same as vb defined in 
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Figure 1. Flux spectral index (/3) distribution of a PWN with 
p = 2.5, nebular radius R = 1 pc, B = 100 fiG, and D = 10 27 cm 2 
s — 1 in a pure diffusion model. 

IGrattonl (Il97l : 

u R = 1 X 10 17 
1 pc 



D 



10 27 cm 2 
100 fiG 



B 



Hz. 



(8) 



We assume a PWN in steady state with p = 2.5, and scale 
to a nebular size i? = 1 pc, magnetic field B — 100 /iG 
and diffusion coefficient D = 10 27 cm 2 s _1 . We calculate 
the spectral index distribution for the three cases: v <C 
vr,v w z/r, and f ^ vr (Figure [1]). When ^ <C vr, 
which corresponds to -q 3> 1, the photon index profile 
is flat. At high frequency the photon index profile first 
changes into a power law when v ss i/r, r\ ss 1, and then 
into an exponential when ^ 3> ^R,f? *C 1. When the 
diffusion distance d 3> R, the diffusing particles within 
the nebula are well mixed so the spectral index profile 
tends to be flat. When the diffusion distance d < R, the 
particle density drops quickly along the radial direction 
because of the short cooling time, so the spectral index 
shows steepening in the radial direction. 

A calculation of the spectral index distribution in the 
KC model shows that it has a problem in explaining 
the o bserved spectral index profile (see also iRevnoldsl 
2003). We take the Crab as an example and use the 
KC model to calculate the spectral index distribution of 
the Crab at optical wavelengths, where the KC model 
is still applicable. We use the best fit parameters given 
bv lKennel fc Coron iti (1984b) to do the calculation and 
assume that there is no synchrotron emission within the 
termination shock. In order to give a good comparison, 
we use the same emissivity as used in lKennel fc Coronitil 
(|1984bl) . which is slighty different from our value. We add 
a pre-defined radius R which is 20 times of the termina- 
tion shock radius in the simulation; there is no emission 
beyond this point. The results are shown in Figure [21 
The spectral index p rofile in Figure [2] does not fit the 
optical data shown in iVeron-Cetty fe; Woltjerl (|1993f ) . In 
the observations, the spectral index profile is approxi- 
mately a power law distribution, while the results given 
by the KC model are flat within a certain radius and then 
increase very quickly beyond that radius. The power law 
like spectral index dis tributions are als o seen at X-ray 
wavelengths in 3C 58 (|Slane et al.ll2004D and G21.5-0.9 
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Figure 2. Flux spectral index (/3) dis tribution of the Crab at opt i- 
cal wavelengths, based on the model of Kennel fc Coronitil i 1984bT ). 
The surface brightness is normalized to the value at the center of 
Crab. 



(jSlane et al.ll2000l; iSafi-Harb et~aT]|2001[ ) which indicates 
that diffusion processes could b e generally impor t ant in 
young PWNe i|Revnoldsl[2003] ). iDel Zanna etall (|2006h 
show that a 2D MHD simulation could reproduce most 
of the toroidal and jet like structure near the termina- 
tion shock, but the spectral index properties of the Crab 
Nebula suggest diffusion processes on larger scales. 

The nebular size of a PWN is also determined by vr 
or r\. When v < Ur, r\ > 1 in the steady state case, the 
nebular size remains the same due to the boundary con- 
dition. When v > vr : r\ < 1, the size tends to shrink as 
the cooling time of particles is smaller than the diffusion 
time. In the v > vr, rj < 1 regime, the nebular size 
for the pure diffusion model can be estimated by setting 
r/ = 1, which yields R ks (6D/QE) 1/2 . 

We first use parameters known from observations such 
as age t, nebular size R and magnetic field B to fit the 
spectral index distribution of the Crab, 3C 58 and G21.5- 
0.9. Fitting a model yields the diffusion coefficient D and 
p of the PWNe. We then discuss the nebular size behav- 
ior. For the Crab we use a magnetic field B = 300 /iG. 
This is slightly l arger than the value, ~ 200 uG, fou nd by 
Ide Jager et all (|1996f l and lAharonian et all <|2004l ) from 
inverse Compton emission, which gives the current value 
of the field. Our value gives a sufficiently high diffu- 
sion coefficient D to explain the size of the Crab. Radio 
data for the Crab show p = 1.52 . The spectral index 
distribution from radio (5 GHz) to optical 6 x 10 14 Hz 
frequencies is shown in Figure [3] By comparin g our re- 
sults with the major axis data in lWilsonl {1972) we find 
that D = 2.5 x 10 26 cm 2 s _1 gives a good fit to the data. 
Si nce we do not know exactly the optical frequency used 
m IWilsonl (fl972h . we did not do a least squares fit. We 
then used the diffusion coefficient D — 2.5 x 10 26 cm 2 s _1 
to calculate the spectral index distribution for infrared 
(3.6 — 4.5 fim) and optical (5364—9241 A) wavelengths, 
which arc shown in Figures [4j In the infrared (IR) , fit- 
ting the spectral index var iation from . 3 to .8 within 
the Crab nebula found by iTemim et alJ (|2006h requires 
the nebular radius to be ~ 130". Our simulation indi- 
cates that in the IR the nebular size of the Crab has 
decreased due t o synchrotron losses , which is consistent 
with Figure 3 of lTemim et alJ (|2006l ) . The spectral index 
variation from 0.6 to 1.1 within the Crab Nebula at op- 
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Figure 3. Flux spectral index (/?) distribution of the Crab from 
5 GHz to 6 X 10 14 Hz, assuming p = 1.52 and B = 300 pG in a 
pure diffusion model. 
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Figure 4. Flux spectral index (/3) distribution of the Crab at IR 
and optical wavelengths, assuming p = 1.52, B = 300 pG, and 



D = 2.5 x 10 26 



in a pure diffusion model, 



tical wavelengths ()Veron-Cettv fc Woltieiifl993l) implies 
the nebular size of the Crab is ~ 100", which is consis- 



tent with the results in Figure 2 of lAmato etall {2000). 
Comparing our simu l ation r esult s in Figure s [3] and HI with 
the da ta in IWilsonl (| 19 721 ) and IVeron-Cettv fe Woftierl 
(1993) shows that our spectral index in the central re- 
gion is lower than the observed value if the optical band 
is involved. We attribute the discrepancy to an intrin- 
sic break in the injected particle spectrum, which is dis- 
cussed below. The main uncertainty in the model comes 
from the magnetic field B, although other factors, such as 
non-spherical symmetry and boundary conditions, also 
have some effect. 

We used the same formalism to fit the X-ray photon in- 
dex profiles of 3C 58 and G21.5-0.9. According to equa- 
tion ([7]), the X-ray emitting particles of both 3C 58 and 
G21.5-0.9 have reached a steady state if they are young 
PWNe, so their age information is not required for pho- 
ton index fitting and we can use the steady state solution 
to calculate the photon index distribution. The parame- 
ters we used for 3C 58 and G21.5-0.9 are listed in Table 
[TJ For 3C 58, the magnetic field B = 80 p,G is based o n 
the minimum energy condition (jGreen fc Scheuerlll99l . 
For G21. 5-0.9, the magnetic fiel d B = 180 /xG is base d 
on the equipartition condition (S afi-Harb et al.l [2001). 
We used a least squar es method to fit the photon in- 
dex data of both 3C 58 (jSlane et al.ll200l and G21.5-0.9 
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Figure 5. Flux spectral index distribution of 3C 58 from 2.2 
keV to 8 keV, assuming p = 2.93, B = 80 [iG and D = 6.11 X 10 27 
cm 2 s _1 in a pure diffusion model. 




Figure 6. Flux spectral index distribution of G21.5-0.9 from 
0.5 keV to 10 keV assuming p = 2.08, B = 180 pG and D = 
3.66 X 10 27 cm 2 s — 1 in a pure diffusion model. 



(|Slane et al.ll2000T ). The 2 parameters in each fit are the 
injected particle spectral index p and diffusion coefficient 
D. The best fit for 3C 58 between 2.2 keV and 8 keV 
with xl e d = °- 83 S ives P = 2 - 93 and D = 6.1 x 10 27 cm 2 
s" 1 (Figure E). For G21.5-0.9, the best fit between 0.5 
keV and 10 keV with X 2 red = 3.28 gives p = 2.08 and 
D = 3.7 x 10 27 cm 2 s" 1 (FigureiJ). The diffusion coeffi- 
cients D we obtained for 3C 58 and G21.5-0.9 are higher 
than for the Crab. Part of the reason may be uncertain- 
ties in the magnetic field B of 3C 58 and G21.5-0.9. As 
discussed above, the diffusion coefficient D and magnetic 
field B follow D cx B 3 / 2 . If the magnetic field B is lower 
than our estimate, the diffusion coefficient D also drops; 
this would imply a high particle energy in the nebulae 
because we used the minimum energy value of B. We 
found that the value of p is insensitive to B and D. 

By using equation (|8|), we obtain i/s = 2x 10 13 Hz for 
the Crab if D = 2.5 x 10 26 cm 2 s" 1 and B = 300 /xG; 
v R = 1.3 x 10 18 Hz for 3C 58 if D = 6.1 x 10 27 cm 2 
s- 1 and B — 80 fiG; v R = 2.6 x 10 17 Hz for G21. 5-0.9 
if D = 3.7 x 10 27 cm 2 s" 1 and B = 180 a*G. For the 
Crab, X-ray, optical and near-IR frequencies are all in 
the v > regime, so the nebular size decreases from 
radio to X-ray. For 3C 58 and G21.5-0.9, all frequen- 
cies below soft X-rays are in the v < vr regime, so the 
radio, optical and soft X-ray nebular sizes of 3C 58 and 
G21.5-0.9 tend to be similar. The different behavior of 
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Table 1 

Parameters used for modeling photon index profile 



Object 


Frequency band 


Magnetic field B 


Distance 


Angular size 


Age 






( M G) 


(kpc) 


of PWN (arcsec) 


(y) 


Crab 


5 X 10 9 - 6 X 10 14 Hz 


300 


2.0 


190 


957 


Crab 


3.6 — 4.5 fim 


300 


2.0 


190 


957 


Crab 


5364 - 9241 A 


300 


2.0 


190 


957 


3C 58 


2.2 - 8 keV 


80 


3.2 


100 




G21. 5-0.9 


0.5-10 keV 


180 


5.0 


40 
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D=2.5x10 26 cm 2 /s, p=1.52 
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Figure 7. Crab Nebula half-light radius based on a pure diffusion 
model with B = 300 fj,G. 



nebular size as a function of frequency among the Crab, 
3C 58 and G21.5-0.9 is due to the fact that the Crab has 
a larger magnetic field but lower diffusion coefficient. We 
use our pure diffusion model to calculate the half light 
radius of the Crab Nebula with p = 1 .52, B = 300 fj,G 
and D = 2.5 x fO 26 cm 2 s^ 1 (Figure [7]), assuming that 
there is an upper limit to the injected particle energy E 
which corresponds to a frequency of 5 x 10 22 Hz. There is 
a bump in the half light radius p lot, mainly becau se for 
the Crab p = 1.52, which is < 2 (jPacholczykll 19701) . and 
we are assuming synchrotron radiation only emits at the 
peak frequency. If we consider the full synchrotron spec- 
trum, the bu mp is diminished. Comparing our results 
to Figure 2 in lAmato et al.1 (|2000h shows that our model 
prediction gives half light radii near the lower limit of ra- 
dio and optical data for the Crab. At X-ray wavelengths, 
our theoretical half light radius is much smaller than the 
observed value. There are several reasons for this. First, 
the spherical symmetry assumption breaks down at X- 
ray wavelengths for the Crab. The Chandra image of the 
Crab s hows clear to roidal structure near the termination 
shock (|Hesterll2008T) . Second, in our model we assume a 
point source while it is in fact an extended source. The 
termination shock has an angular size ~ 10" (0.1 pc) at 
the Crab Nebula. At radio and optical wavelengths, the 
nebular size is much larger than the size of the injection 
region so the point source assumption is adequate, but at 
X-ray wavelengths it is no longer true. The last reason is 
that our pure diffusion model does not include the effect 
of advection. It is likely that both advection and diffu- 
sion play a role in PWNe. Assuming V a d v oc r~ 2 near the 
pulsar wind termination shock (KC), t a d v /tdiff oc R, so 
advection becomes more important in the inner regions. 
We expect advection to play some role in X-ray emission 
from the Crab, and we discuss it in Section 3.2. 
In considering the integrated number density N(E,t), 



we no te the result for synchrotron losses only ([Pacholczvkl 
1970) 



N(E,t) = 
K 

(P-1)Q 
K 

(p-i)Q 



£-(p+i)[i _ (i_ QEty- 1 ] 

£H p+1) , \iQEt>l. 



if QEt < 1 



(9) 



Our pure diffusion model deviates from this result be- 
cause we have a pre-defined PWN radius R and assume 
a transparent outer boundary at R; particles that diffuse 
out of R are not taken into account in the integrated 
number density N(E,t). The results shown in equation 
([9]) would apply to a pure diffusion model with a constant 
magnetic field B and a reflecting outer boundary which 
counts all the injection particles. As discussed in Section 
2, the issue of transmission through the outer boundary 
of a PWN is uncertain from the observational point of 
view. We discuss the spectral index distribution and half 
light radius of a model with a reflecting outer boundary 
in Section 3.2. 

The Crab, 3C 58 and G21.5-0.9 show flat radio spec- 
tra and cannot be fitted by a single power law injec- 
tion spectrum at radio through X-ray wavelengths even 
taking into account the ev ol ution of the PWN (e.g. , 
iRevnolds fc Chevalierl ll~98l . iBucciantini et al.l (|2011[ ) 
considered a 1 zone model with a broken power law in- 
jection spectrum and long term evolution of the PWN, 
showing that it can explain the integrated spectra of the 
Crab and 3C 58 from radio to X-rays. The intrinsic spec- 
tral b reaks for all PWNe c o nside red are at a similar en- 
ergy. iSironi fc Spitkovskvl (|2011[ ) did both 2D and 3D 
particle in cell simulation for the termination shock and 
show that it could create both a flat power law (p ~ 1.5) 
and steep power law {p ~ 2.5) components in the post- 
shock spectrum. The spectral break of PWNe between 
radio and optical wavelengths may be a natural conse- 
quence of particle acceleration at the termination shock. 
In fitting the optical spectral index distribution for the 
Crab, we already mentioned that there are indications of 
another power law component in the Crab. Here we use a 
double power law injection spectrum in our pure diffusion 
model and re-calculate the spectral index distribution 
and half light radius of the Crab. The evolution of the 
PWN is still ignored here because of the additional com- 
plications. For 3C 58 and G21.5-0.9, the X-ray spectral 
index data are well above the spectral break frequency 
and our modeling is not affected by the double power law 
feature. We continue using a magnetic field B = 300 //G 
and p\ = 1.5 for the low energy component of Crab par- 
ticle distribution. We use P2 = 2.35 and a break energy 
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Figure 8. Flux spectral index (/3) distribution of the Crab for 
different wavelength bands assuming pi = 1.52, p2 = 2.35, B = 300 
fiG, and D = 6.0 X 10 26 cm 2 s — 1 in a pure diffusion model. 

4 x 10 11 eV for the o t her p ower law component as found 
bv lBucciantini et all (|2011D and find that D = 6. x 10 26 
cm 2 s " 1 gives a good fit to the major axis data in I Wilson] 
(|1972f) , as shown in Figure [5] Then we use the same dif- 
fusion coefficient D to calculate the spectral index distri- 
bution at IR and optical wavelengths (Figure [5]) ■ After 
adding another power law component, w e obtain a bet- 
ter fit to the central spectral index data in I Wilson] (|1972f ) 
and iVeron-Cettv fc Woltierl (|1993l ). The nebular radius 
we need to explain the spectra l index variation at IR 
(iTemim et alj |2006.) and optical (|Veron-Cettv fc Woltierl 
[1991 " wavelengths becomes smaller: ~ 110" in IR and 
~ 70" in optical wavelengths. Our simulation gives a 
spectral index of 0.7 at the center of the Crab at optical 
wavelengths, which is slightly larger than the observed 
value. This is due to t he fact that we use p 2 = 2.35 
(jBucciantini et aH 1201 If ) fo r the steep power l aw. H ow- 
ever, the best fit parameter iBucciantini et"aT1 (|2011| ) ob- 
tained in their ID evolution model may not be the best 
fit parameter for our case as we are considering a steady- 
state situation. If we change the break energy for the 
two power laws, we could obtain a better fit to the opti- 
cal nebular radius, but the improvement is not significant 
in view of the other uncertainties in the model. The half 
light radius of the Crab with double power law injection 
spectrum is shown in Figure [7] The double power law 
fit gives a larger nebular size in the radio band which is 
in the non-steady state regime, mainly because we use a 
larger diffusion coefficient D in the double power law fit. 

So far, the results are under the assumption that all 
the emission of an electron goes into radiation at a fre- 
quency v corresponding to the maximum synchrotron ra- 
diation power. We will continue to make this assumption 
in the next section because it speeds the calculations, but 
here we carry out the calculation with a full synchrotron 
spectrum to show the uncertainty caused by our approx- 
imation. We considered the Crab, 3C 58 and G21.5- 
0.9, and our numerical results with the full synchrotron 
spectrum (Figure [5]) show that the diffusion coefficient D 
required to fit the observations drops by a factor about 
2. Here, we in tegrate the synchr otron radiation function 
F(x = v/v c ) ()Pacholczvklll970D from 0.005 to 500 and 
assume that for 3C 58 and G21.5-0.9 all the particles 
within that energy range are in a steady state. Because 
F(x = v/v c ) drops very quickly beyond the peak and 



we are calculating the spectral index distribution of 3C 
58 and G21.5-0.9 in X-rays, time dependent effects are 
expected to be small. 

In the pure diffusion model we have assumed that the 
diffusion coefficient is constant. However, it is possi- 
ble that the diffusion coefficient in the PWN has energy 
dependence De(E), as discussed in Section 2, and we 
now investigate how the energy dependence of the dif- 
fusion coefficient would affect our pure diffusion model 
with transmitting boundary. We use the Green's func- 
tion method to solve the steady state equation of the 
pure diffusion model but now with an energy dependent 
diffusion coefficient De{E) = D E a : 



D E a V 2 N + Q 



dE 2 N 
dE 



-KE- p S(r) 



(10) 



The resulting particle distribution function is (see Ap- 
pendix) 



N(E, r) 



K 



4:TrrD E a ' 



-E-P- 



1 



where 



(l ) 1 —=dx, 

V x) Jx 



2 QE{l-a) 



(11) 



AD a E a 



In the solution, u changes to u — r 2 QE(l — a)/4D E a ; 
again, ^/u can be considered as a ratio of the diffusion dis- 
tance to the nebular radius. In a steady state, t = 1/QE, 
and we have y/u = r(l - a) 1/2 /(4D S Q t) 1/2 . Setting 
u = 1, we find that the nebular size R cx E 1 ^ 1- ")/ 2 oc 
i^ 1- ")/ 4 . For the spectral index distribution, we note 
that a > implies that more energetic particles dif- 
fuse out more rapidly than less energetic particles, which 
should flatten the spectral index distribution; the same 
effect results from reducing the magnitude of the diffu- 
sion coefficient when a — 0. Our simulations show that 
the case De = DqE" corresponds to a constant diffu- 
sion coefficient case with D if De = (1 — a) 2 D, where 
De is taken at the energy at which the spectral index is 
measured. The flux spectral index at the center now is 
(pe + a — l)/2 (equation [11]) because the integral part 
in the solution approaches some constant when r — > 0. 
In order to get the same spectral index at the center as in 
constant diffusion coefficient case we need pe — p— a. In 
Figure HO] we consider a PWN withps = p — a = 2.5 — a, 
nebular size R = 1 pc, magnetic field B = 100 /iG, and 
diffusion coefficient De — 10 27 (1 — a) 2 cm 2 s _1 at an 
energy corresponding to v — 10 17 Hz, and plot the spec- 
tral index profile for different a values. It is clear that 
they all have a similar slope, implying that for a cer- 
tain band, if a pure diffusion model with constant D can 
fit the data, a pure diffusion model with energy depen- 
dent De can also fit the data but with different p and 
diffusion coefficient at that band. Although energy de- 
pendence of the diffusion coefficient does not change the 
consequences of spectral index fitting with constant D at 
one wavelength band, it changes the fitting result when 
multiband data are considered because now the nebular 



size R oc v a '/ 4 instead of R 



oc v 



-1/4 



We previously 



found with the constant diffusion coefficient assumption, 
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Figure 9. (a) Flux spectral index (/?) distribution of Crab from 5 
(b) Flux spectral index (fi) distribution of 3C 58 from 2.2 keV to 8 
distribution of G21.5-0.9 from 0.5 keV to 10 keV assuming p = 2.08 

there was a problem in explaining the nebular size of the 
Crab beyond optical frequencies. If we take the energy 
dependence of Dg into account, we can fit the data on 
nebular size and overall appearance of the spectral in- 
dex profile in the IR and optical better. However it is 
difficult to determine whether energy dependence is re- 
quired by the data. Since other factors like advection 
affect the nebular size and spectral index profile of the 
Crab in a similar way, and the observational data for 
the Crab Nebula do not have high precision, the energy 
dependence of the diffusion coefficient is not determined 
by our models. However, the decreasing size at higher 
photon energies implies that a < 1; as discussed in Sec- 
tion 2, some treatments of diffusion in old PWNe have 
assumed a = 1 . 

3.2. Diffusion and advection model 

We now use Monte Carlo simulations to consider mod- 
els with diffusion and advection. Monte Carlo methods 
allow the treatment of cases that cannot be treated in 
the pure diffusion, analytical model. We expect that 
advection plays some role in the spectral index distri- 
bution and half light radius, especially in the v vr 
regime. In the MHD model of KC, there is an advective 
flow after the termination shock in which the flow ve- 
locity declines from mildly relativistic to the velocity at 
the edge of the nebula, ~ 2000 km s _1 . If the magnetic 



GHz to 6 X 10 14 Hz assuming p 1 = 1.52, p 2 = 2.35 and B = 300 fjtG, 
keV assuming p = 2.93 and B = 80 /iG, (c) Flux spectral index (/}) 
and B = 180 (jG. 
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Figure 10. Flux spectral index (/3) distribution at 10 17 Hz for 
a PWN with pe = 2.5 — a, nebular size R = 1 pc, magnetic 
field B = 100 /LtG, and diffusion coefficient D E (v = 10 17 Hz) = 
10 27 (1 — a) 2 cm 2 s -1 , for various values of a. A steady state 
model with transmitting boundary is assumed. 

field is not dynamically important, the velocity declines 
as r~ 2 in a steady flow. Another consideration is a re- 
flecting outer boundary condition. For PWNe satisfying 
r\ = (6Dt/ 'i? 2 ) 1 / 2 3> 1, the reflecting boundary should 
have a significant effect on the spectral index distribu- 
tion and half light radius. For the reflecting boundary 
effect we primarily focus on 3C 58, which has a high dif- 
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fusion coefficient D. The last improvement is treating 
an extended source, since the termination shock has a 
finite size ~ 0.1 pc which can play a role in the cen- 
tral region. The time dependence of the PWN and non- 
s pherical symmet ry are not taken into account. 

iMassarol (|1985l ) discussed a diffusion and advection 
model, but many assum ptions were needed to ob- 
tain an analytical mode l. IWcinbcrg fc Silkl (|1976l ) and 
iRevnolds fc Jones] (|1991|) investigated approximate solu- 
tions for a pure diffusion model with an extended source, 
but they omit some of the physical effects that we wish 
to include. We therefore developed a code to carry out 
Monte Carlo simulations that allowed a general treat- 
ment of diffusion and advection processes. We modeled 
both the Crab and 3C 58. For 3C 58 we mainly focus 
on the variation of spectral index distribution while for 
the Crab we are interested in the variation of the half 
light radius with photon energy. In simulating 3C 58, a 
total number of 10 6 effective particles were injected into 
a spherical shell at time intervals of 5.4 x 10 4 s. Since we 
are considering the X-ray band of 3C 58, which is in a 
steady state, only particles injected within a time 1/QE 
of the present are taken into account. In simulating the 
Crab, at least 10 6 effective particles were injected into a 
spherical shell at time intervals of a half day. The motion 
of the particles is a superposition of advective motion and 
3-dimensional random motion. The displacement in each 
time interval is 

AStct = V adv Atr ± (2D x At)^ 2 x Q 
±{2D y Mf/ 2 y ± (2D z At) 1 ' 2 z a , (12) 

where x , y a , z and r are unit vectors in x, y, z and ra- 
dial directio ns. We take V a ( j v oc r~ 2 , V(r = R out ) = 
630 km s~ 1 (|Bietenh olz 2006j) and the ratio of outer ra- 
dius to in ner radius to be R ou t/Rin = 100/8 = 12.5 
for 3C 58 (ISlane et al.H2004D . In this case, we only con- 
sider the synchrotron radiation loss since the adiabatic 
expansion energy loss, E oc V • V a dv, is zero for our 
velocity profile. For the Cr ab, we use some of the pa - 
rameters from the model of IKennel fc Coronitil ([T984b), 
assuming that the velocity at the termination shock is 
c/3, the velocity V a dv oc r~ 2 between the termination 
shock radius R s and a critical radius R c , and V a dv re- 
mains constant between the critical radius R c and the 
outer radius of the ne bula R ou t- We ass ume i? s = 15 
arcsec in angular size (jHester et al.ll200"2l) . and in order 
to set V a dv (Rout) = 2000 km s _1 (jKennel fc Coronitil 
Il984aft we further assume R c = 50 1 / 2 R S . We consider 
adiabatic expansion energy losses in addition to syn- 
chrotron radiation losses for R c < R < R ou t\ it is 
dE/dt = —2V a dv(Rout)E/3r for relativistic particles in 
a medium with constant flow velocity. In all the simula- 
tions, we use a constant magnetic field B to simplify the 
calculation. In spherical symmetry, D x = D y = D z = D. 
If particles move into the inner boundary due to their 
random motions, they are forced to bounce back into 
the PWN region. A reflecting boundary was used at the 
outer radius. 

The simulation results for 3C 58 are shown in Figure 
ITT1 The additional physical processes allow a better fit 
to the data. In order to discern what role advection 
and a reflecting boundary play in the photon index dis- 



tribution, we did one simulation for pure diffusion with 
only an inner reflecting boundary and one for pure dif- 
fusion with both inner and outer reflecting boundaries 
(Figure ITT]) . After comparing the results of the two cases 
(Figure [UJ, we find that the outer reflecting boundary 
condition mainly makes the photon index steeper and 
the part near the outer boundary relatively flat. The 
diffusion coefficient D required to fit the data becomes 
higher. Advection is not very important in the fit be- 
cause the ratio of diffusion time to advection time ratio 
is low. In the Figures II H we have shown the ratio of 
diffusion time to advection time tdi / / /t a d v , where t d i / / 
is estimated by (61?^//) 1/2 = R ou t ~ Rin and t adv is 
calculated by t a d v = J^."" 4 dr/V a dv The angular size of 
the flat region in the radial direction can be estimated as 
follows. Since tdiff <C t a d v in our simulation, advection 
can be neglected. For particles in a steady state, the dif- 
fusion distance is R dlff = (QDt) 1 / 2 = (6D/QE) 1 / 2 , so 
that 9 fi a t ~ 6 x (Rdiff/R — 1), where 9 is the angular 
size of PWN. Substituting the expression for Rdiff and 
recalling that En satisfies R = (AD/QEr) 1 ' 2 , we obtain 
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For 3C 58 with D = 8.8 x 10 27 cm 2 s" 1 , B = 80 fiG, and 
R = 1.55 pc, equation (JT3J) implies 8fi at ~ 30" which 
is consistent with the results in Figure [UJ The results 
shown in Figure [TT] are calculated for 2.2 — 8 keV pho- 
ton energies. When we do the calculation for 9fi a t we 
use 8 keV as it gives a lower diffusion distance. We em- 
phasize that equation (fT3|) is only correct for a steady 

state and requires that 1 > ^(^)^ -1 > 0. When 

2 C^)* — 1 < 0, the cooling time is lower than the dif- 
fusion time, and few particles reach the outer boundary. 

When - 1 > 1, due to the boundary effect, the 

diffusion coefficient D is large enough to smear out the 
spectral index structure in the nebula, so the spectral 
index distribution is flat within that energy band. 

We calculate the spectral index distribution of the 
Crab at optical wavelengths (Figure [T2]) and its half light 
radius from ultraviolet (UV) to X-ray frequencies (Ta- 
ble [H) , which are above the break frequency given by 
iBucciantini eTaTl (f20lTh . with p 2 = 2.0, B = 300 fiG 
and D = 9.0 x 10 26 cm 2 s _1 . The spectral index vari- 
ation from 0. 6 to 1.1 within the Crab Neb ula at optical 
wavelengths (jVeron-Cetty fc Woltierl["l993! ) now implies 
the nebular size of Crab is ~ 140", which is slightly larger 
than the pure diffusion case as we now use a larger diffu- 
sion coefficient D and take advection into account. The 
advection process increases the half-light radius of the 
Crab significantly at X-ray energies. Based on Table [2j 
the ratio (Rdif / - Rs)/ (Rdiff +adv~ Rs) increases rapidly 
from X-ray to UV wavelengths as we assume the ex- 
tended source size is ~ 15". The reason advection is im- 
portant for the Crab in X-rays is because the X-ray size of 
the Crab is small. As noted in Section 3.1, the ratio of ad- 
vection time to diffusion time R a dv/ Rdif f oc R, the neb- 
ular size; advection becomes more important when the 
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Figure 11. (a) Monte Carlo simulation with diffusion and advection for 3C 58 from 2.2 keV to 8 keV assuming p = 2.8, B = 80fiG, and 
reflecting inner and outer boundaries, (b) Monte Carlo simulation with pure diffusion, a reflecting inner boundary, and a tranmitting outer 
boundary, (c) Monte Carlo simulation with pure diffusion and reflecting inner and outer boundaries. 

1.5 | — . — , — , — , — , — , — , — , — , — , — , — , — , — , — , — , — , — , — 

Table 2 

Half light radius of the Crab Nebula in arcsec 



Frequency 


10 15 Hz 10 16 Hz 


10 17 Hz 


10 18 Hz 


pure diffusion 


72 43 


27 


19 


diffusion and advection 


84 55 


35 


24 



nebula is small. The half-light radius derived from the 
diffusion and advection model can now ex plain the high 
frequ ency part of the Crab nebula size data ([Amato et al.l 
2000). However, the half-light radius we obtained still 
drops sharply as a function of frequency, which may im- 
ply energy dependent diffusion. The values of P2, mag- 
netic field B, and diffusion coefficient D we choose are 
not best fit parameters. Certain combinations of parame- 
ters would improve the fit to the observational data. The 
velocity and magnetic field profiles for the diffusion and 
advection case must be analyzed for more exact models. 
The uncertainty in the termination shock radius and the 
flow velocity at the outer radius of the nebula also affect 
the simulation results. 

As discussed in Section 1, 2-dimensional MHD mod- 
els reproduce many features observed in the inner Crab 
Nebula. Diffusion is not a factor in this region for 2 rea- 
sons: the short advection timescale because of the high 
flow velocities and a long timescale for radial diffusion 
because of the toroidal magnetic field. In the Crab, the 
prominent toroidal structure observed at X-ray and opti- 



5364 A to 9241 A 



0.5 



100 

radius(arcsec) 



Figure 12. Flux spectral index (J3) distribution of the Crab at 
optical wavelengths assuming p = 2.0, B = 300 fiG, and D = 
9 X 10 2e cm 2 s in a Monte Carlo simulation with reflecting inner 
and outer boundaries. 

cal wavelengths extends to ~ 40" from the pulsar, while 
the nebular radius is 200". Our model applies to the 
outer 4/5 of the nebula. Toroidal structure is less promi- 
nent in 3C 58 and G21. 5-0.9. 

We also used Monte Carlo simulations to investigate 
the case where there is particle transport across the neb- 
ular boundary at R and particles are lost from the system 
once they cross R. The effect is to flatten the spectral 
index profile in the outer part of the nebula and to re- 
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duce the value of D by a factor of 2 compared to the 
simple model (Section 3.1). However, we consider the 
reflecting boundary model to be more plausible, for the 
reasons given in Section 2. 

4. DISCUSSION AND CONCLUSIONS 

We have argued that the structure of young PWNe is 
not described by a tor oidal magnetic field , as exp ected 
in a model like that of Ke nnel fc CoronitH (fi~984a), but 
has a more chaotic magnetic structure. In Section 2, 
we noted various observational studies that showed con- 
siderable structure in young nebulae, not a clear toroidal 
structure. In Section 3, models with diffusion of particles 
were presented and compared to observations of 3 neb- 
ulae. Emphasis was placed on fitting the spectral index 
profiles of the nebulae, as well as the surface brightness 
profiles. Models with diffusion were much better able to 
spectral index profiles than pure advection models. The 
best estimates of the diffusion coefficient come from the 
Monte Carlo simulations, but these values need to be 
somewhat reduced because they do not include the full 
synchrotron spectrum in the calculation of the emission. 
Estimates of the diffusion coefficient and the correspond- 
ing particle mean free path are given in Table [H The 
assumed magnetic field is given for each case because 
D cx B 3 / 2 and the magnetic field strength is uncertain. 

Table [3] shows that the diffusion coefficient and length 
for the Crab is considerably smaller than that for 3C 58 
and G21.5-0.9. The length does not scale with the size 
of the PWN because the Crab is larger than both 3C 58 
and G21.5-0.9. One possibility is that there is frequency 
dependent diffusion coefficient. The coefficient for the 
Crab is derived from optical/IR observations, so that the 
lower energy particles are being observed compared to 
3C 58 and G21.5-0.9 where X-ray observations are used. 
Table |3] shows the corresponding particle energies for dif- 
fusion coefficients. An energy dependent diffusion with 
a ~ 0.5 — 0.6 would explain the difference between the 
Crab and the other remnants. As discussed in Section 
3.1, this is consistent with our results and we suggest it 
as a possible explanation for the difference between the 
PWNe. 



Table 3 

Diffusion coefficient and length 



Object 


D 


B 


A = D/c 


Particle energy 




(cm 2 /s) 


0*6) 


(10 16 cm) 


(TeV) 


Crab 


2.4 X 10 26 


300 


0.8 


0.6 


3C 58 


2.9 X 10 27 


80 


10 


40 


G21.5-0.9 


2.0 x 10 27 


180 


7 


30 



The spectral indices along magnetic filaments in the 
Crab shows relatively little variation while the spectral 
indices show steepe n ing away fro m the filament center 
(|Seward et all 120061: iHesterl 120081) . which is consistent 
with rapid particle motions along filaments and slow dif- 
fusion across filaments. However, the diffusion mean free 
paths (Table are smaller than characteristic structures 
in the nebulae. The length for the Crab is about a factor 
10 smaller than the scale indicated by optical polariza- 
tion (Section 2), and the length for 3C 58 is about a 
factor of 10 smaller than the sc ale of apparent ma gnetic 
loops seen in the X-ray image (jSlane et al.l I2004D . The 



actual longer diffusion time (due to the smaller length) 
may indicate that the magnetic structure is not com- 
pletely random. 

Our models have been designed for comparison with 
young PWNe like the Crab, 3C 58, and G21.5-0.9. 
They are likely to be in an evolutionary phase where 
the nebulae are accelerating into freely expanding su- 
pernova ejecta. In a subsequent phase of evolution, 
the reverse shock wave from interaction with the in- 
terstellar medium comes back toward the center and 
can push off the PWN, creating an asymm e tric n ebula 
(jBlondin et al.|[200l . IVan Etten fc Romanil (pOll have 
investigated evolutionary models for the PWN HESS 
J1825-137, which probably belongs to the class of post- 
reverse shock nebulae and is observed at X-ray and TeV 
energies. Of interest is the fact that their modeling 
shows the ne ed for diffusion of p articles. As mentioned 
in Section 2, IHinton et al.l (| 2 1 If) have considered diffu- 
sion from the evolved PWN Vela X. A chaotic magnetic 
field is expected in these objects because of instabilities 
related to the interaction wit h the reverse shock front 
from the supernova remnants (jBlondin et al.| [2001). The 
situation is different for the young remnants discussed 
here, but a chaotic field may be the res ult of instabilities 
related to the toroidal magnetic field ()Begehnanl 119981) 
and Rayleigh- Taylor instabilities in the outer parts of 
the nebulae. 
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APPENDIX 
In order to solve the equation 

D E a V 2 N + Q-g^- = -KE~P6(f), (14) 

we first let $ = QE 2 N and 4nJ = ^E 2 - a KE-PS(r), 
so equation can be simplified to 



Q 



d<S> 



V^--g-(l-a)^ ZT = -4 7 rJ. (15) 



dE 



a-l 



WO 



(16) 



Further, assuming m 2 = -^(1 — a) and r = E 
have 

9 9 9$ 

or 

Here we only consider the case m 2 > 0, which requires 
a < 1, as it is consistent with the situation in the Crab. 
If a > 1, the nebular size of the Crab does not decrease 
with increasing frequency. Next we consider the Green's 
function equation corresponding to the above equation: 



V 2 G- m 2 



dG 
dr 



-4ir8(T-T')8(r-P). (17) 



If we do not consider the boundary effect, the Green's 
function solution for equation (|17|) is 



G(r, T /r',T') 



rr> 



2^(t-t') 1 ' 5 



Ht-t') u (j — r ') j 

(18) 



13 



where u(t — r') is the step function. Here, G(r, r /r 1 ,r') 
is already normalized to 

777^ /* -> 

Then 

$( r ,r) = ^ dr' J dV'J{r',T')G(r,T/r',T'). (20) 
After some calculation, we obtain 

v ' ; D 47r 1 - 5 r J u \ xJ y/x ' 

(21) 

where w = Q(1 4 ~° )r2 ff 1 -". Since iV(r,£7) = <S>/QE 2 , we 
finally obtain 

(22) 
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